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We reformulate the strong-interaction limit of electronic density functional theory in terms of 
a classical problem with a degenerate minimum. This allows us to clarify many aspects of this 
limit, and to write a general solution, which is explicitly calculated for spherical densities. We then 
compare our results with previous approximate solutions and discuss the implications for density 
functional theory. 
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I. INTRODUCTION 

Density functional theory [1-3] (DFT) is by now the 
most popular method for electronic structure calculations 
in condensed matter physics and quantum chemistry, be- 
cause of its unique combination of low computational cost 
and reasonable accuracy for many molecules and solids. 

In applying DFT to a given electron system, the only 
quantity that must be approximated in practice is the 
functional -E X c[p] for the exchange-correlation energy. 
An exact expression for this functional is the coupling- 
constant integral, 



The integrand is defined as 



daW a [p}. (1) 



W a \p] 



\V ee \* a [p])-U[p}. (2) 



Here, U[p] = § / dr J dr' p(r)p(r')/\r - r'| is the func- 
tional of the Hartree energy and the operator V ee de- 
scribes the Coulomb repulsion between the N electrons, 



V RR = 



N-l N 

E E ■ 



i 



(3) 



(Atomic units are used throughout this work.) Even- 
tually, out of all antisymmetric TV-electron wave func- 
tions ^> that are associated with the same given electron 
density p = p(r), *& a [p] denotes the one that yields the 
minimum expectation of the operator T + aV C0 . Here, 
T = -^Eili i s the kinetic-energy operator. Notice 
that the parameter a works as an adjustable interaction 
strength or coupling "constant" . 

At a — 0, W a [/o] starts out with the value 



E x [p] = (* [p]\V ee \* [p})-U[p] 



(4) 



which is the density functional for the exchange energy. 
Generally, W a [p] is a monotonically decreasing function 
of a, since the electrons in the state ^a[p] increasingly 



tend to avoid each other in space as the repulsion strength 
a grows. Thus, the expectation of V eei which is a measure 
for the average inverse distances between the electrons, 
must decrease. However, it cannot decrease indefinitely, 
since the electrons are forced to stay within the fixed 
density p(r). 

Schematically, the traditional quantum chemistry ap- 
proach to electron correlation often consists in trying 
to extrapolate the information on the physical system 
(q = 1) from the non-interacting limit (a = 0) by using, 
e.g., perturbation theory or more sophisticated methods. 
This work, instead, follows the early idea of Wigner [4], 
further developed in the DFT framework in Refs. 5 and 
6, in which the information at a = 1 is obtained by in- 
terpolating between the two limits of weak interaction, 
a — > 0, and infinitely strong interaction, a — > oo. In this 
context, we carry out a detailed study of the limit 



Woolp] = lim W a [p] 



(5) 



Although the limit of Eq. (5) was investigated in previous 
work [5, 7-9], the solution presented here was found only 
in the special case of two electrons in a spherical density, 
using physical arguments. For the general case with TV 
electrons, the point-charge-plus-continuum (PC) model 
was proposed [5, 8], 



W^[p] 



dr 



Ap(r) 4 ^ 3 



D 



|Vp(r)| 
p(r)4/3 



(6) 



where A = _^(*p)V3 and B = -^{±yr\ This ap 
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proximation, together with a similar one for the next 
leading term, W^Mj describing zero-point motion oscil- 
lations, was used to construct an interpolation for W a [/>] 
between a = and a = oo, called interaction strength 
interpolation (ISI) model [6]. ISI predicts accurate at- 
omization energies (with a mean absolute error of 3.4 
kcal/mole) [6, 10], showing that the general idea of in- 
terpolating between the weak- and the strong-interaction 
limits in DFT can work. 

It should be emphasized that ISI uses as ingredients 
exclusively the two functionals E x [p] and E^ L2 [p] (of 
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the second-order correlation energy in Gorling-Levy per- 
turbation theory [11]) from the relatively simple non- 
interacting limit a — > (with single-particle orbitals), 
plus the two functionals [p] and [p] from the op- 
posite a —* oo limit of infinitely strong repulsion, which 
are, so far, not known exactly (as said, except for WooM 
in the special case of two electrons in a spherical density 
[7])- 

In this work we reformulate the a —* oo limit of DFT 
in terms of a classical problem with a degenerate mini- 
mum. This allows us to construct the general solution 
for Woo[p], and to clarify many aspects of this limit. As 
a > grows beyond its real- world value a = 1, the con- 
cept of single-particle orbitals becomes completely mean- 
ingless. When a — > oo, however, an entirely new type 
of simplicity with so-called "co-motion" functions arises, 
as we shall see in Sees. II-IV. The "co-motion" func- 
tions, which entirely determine Woo[p], can be directly 
constructed from the one-electron density p(r). 

The paper is organized as follows. In Sec. II we first 
give a general overview of the problem, anticipating the 
solution on intuitive physical grounds, and leaving the 
mathematical details in Sees. Ill and IV. We then use 
our formalism to explicitly evaluate the limit of Eq. (5) in 
the case of spherical densities, with applications to atoms 
(Sec. V). In Sec. VI wc compare our solution with the ap- 
proximation of Eq. (6), and we discuss the implications 
for the ISI functional. The last Sec. VII is devoted to 
conclusions and perspectives. In the Appendix we also 
consider the simple case of harmonic forces, in order to 
analyze how the nature of the electron-electron interac- 
tion affects the solution. 



II. SMOOTH DENSITIES FROM A CLASSICAL 
PROBLEM 

For a given V-electron density p = p(r) we generally 
wish to find the limit (5) or, equivalently, 



W^p] + U[p] = lim (V a {p]\V ee \y c 

a — >-oo 



M>- (7) 



If the density p is both N- and w-representable for ev- 
ery a, there exists an a-dependent external potential 



KxtW = Eili «"xt(bL r *) such that is the ground 



state of the Hamiltonian 



H a [p]=f + aV cc + V^ t [p]. 



(8) 



As a — > oo, the binding external potential v" xt ([p],r) in 
the Hamiltonian (8) must compensate the strong repul- 
sion 0(a) between the electrons. Therefore, we expect 



a — *oo Ol 



(9) 



with v([p], r) a smooth function of r. Thus, for large 
a » 1, the kinetic energy in the state ^f a [p] is mainly 
due to zero-point oscillations of strongly repulsive elec- 
trons, bound by a strong attractive force that has the 



order of 0(a). Therefore, (* a [p]|T|* Q [p]) has the order 
of O(y/o) and we may write in Eq. (7) 



WovW + Ulp] = lim (* a [p] 



1 



-T + V C( 



lim min ( ^ 

a— -co ty^p 



T + V a 



*a[p] / 

*y (io) 



In the second step we have applied the definition of the 
wave function ^ a [p]; the constraint — > p" addresses 
all those wave functions W that are associated with the 
same given density p. Provided that the limit a — > oo 
can be applied directly to the operators in Eq. (10) (there 
is no rigorous proof for this reasonable conjecture - see 
also the Appendix), it is simplified to 



W^ + Uip] 



min (1- ivy*). 



(11) 



In this case, the expectation of V cc alone is to be mini- 
mized, regardless of the kinetic-energy operator T. This 
apparently purely classical problem corresponds to the 
quantum- mechanical limit of infinitely large masses. 

If Eq. (9) holds, that is if the density p(r) is both 
N- and w-representable also in the very a — > oo limit, 
the minimizing * in Eq. (11) is the ground state of the 
pure multiplicative operator V cc + V. The local one-body 
potential V is the Lagrange multiplier for the constraint 
"\f r — > p" , and can be found via the Legendre-transform 
formulation of Eq. (11) [12], 



vMp] + c/[ p ] 



max ^ min(\l/|I4 e 



p v 



(12) 

To start to address the problem of V-representability 
in the a — > oo limit, wc write the minimizing wavefunc- 
tion * in Eq. (12) as the product 



* = -0(ri, ...,rjv)x(<7i, ...,<tn), 



(13) 



where tp is a spatial wavefunction and x is a chosen eigen- 
state of the total spin S 2 and its projection S z . In Sec. Ill 
we analyze the V-representability problem in more de- 
tail, showing that, in the special a — > oo case, we can al- 
ways construct an antisymmetric * of the form (13), and 
that the choice of the spin cigenfunction \ does not af- 
fect the square of the spatial wavefunction \tp{ri, r 2 )| 2 
and thus, as shown by Eq. (15) below, the energy. For 
this reason, in what follows we only consider the spatial 
wavefunction ip(r\, rjv)- 
As said, 



V cc + V ee E pot ([v];r ll ...,r N ) 



(14) 



is a pure multiplicative operator, so that, after integrat- 
ing out the spin variables, the unconstrained minimiza- 
tion in the brackets of Eq. (12) reads 

min / dr 1 ...dr N \ip(r ll ...,r N )\ 2 E pot ([v];r ll ...,r N ). (15) 

■>l> J 
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We see that the minimum in Eq. (15) is reached when 
the square of the spatial wavefunction \tp{ri, ...,i"at)| 2 is 
a distribution that is zero everywhere except for values 
(ri, rjy) G M, where 




, -,rjv) : E pot ([v];r-i, ...,r N ) = min [ (16) 



is the set of all configurations (ri, ...,rjv) for which the 
37V-dimensional function E pot of Eq. (14) assumes its 
absolute (global) minimum. Notice that the set M is 
purely determined by the choice of v(r), M = M[v). We 
shall see in the next Sec. Ill that to such a minimizing 
distribution \ip(ri, r N)\ 2 we can indeed associate an 
antisymmetric wavefunction of the form (13). 

For a given reasonable (attractive) potential v(r), we 
may expect that M will comprise only one single config- 
uration (plus its permutations) or a small finite number 
of configurations: the set M corresponds, in fact, to the 
solution of the classical electrostatic equilibrium prob- 
lem for N identical point charges in the external poten- 
tial v(t). The density associated with the corresponding 
spatial wave function tp will then be a finite sum of delta 
functions, 

K 

p(r)cx^(5(r-r fe ), (17) 

fe=i 

where the are the K > N different position vec- 
tors forming the configurations in M. When u(r) has 
some symmetry (e.g., spherical), M can also be a low- 
dimensional continuum (e.g., the surface of a sphere). 

For a quantum-mechanical system (atom, molecule, 
solid), the electronic density is typically smooth, so that 
Eq. (17) seems in contradiction with the final require- 
ment "\& — > p" [or the maximization (12)]. In other 
words, it seems like there is a v-representability prob- 
lem for smooth densities in the very a — > oo limit. The 
solution to this apparent problem is that the potential 
v(r) must be such that M[v] is a continuum with at least 
three dimensions. That is, the absolute minimum of the 
3iV-dimensional function E pot must be degenerate over 
(at least) a 3-dimensional subspace of R 3N , which we 
write as 

M= {(r,f 2 (r),...,f w (r)) : re?}, (18) 

where P C R 3 is the region of space where p(r) ^ 0. As 
we shall see in Sec. IV, the ansatz (18) can be consistent 
with the requirement that E pot is minimum over M only 
when the functions f„(r) satisfy special properties. 

From the physical point of view, the corresponding spa- 
tial wavefunction ip (which is zero everywhere except on 
M), describes a state in which the position of one of the 
electrons can be freely chosen in P, ri = r, but it then 
fixes the positions of all the other electrons through the 
functions f n (r), r 2 = f 2 (r ),..., rjv = fjv(r). This is what 
we call [5, 7] "strictly correlated electrons" (SCE). 



The requirement "\& — > p" can of course be fulfilled 
also if the degeneracy of the absolute minimum of E pot 
is higher, e.g., if M is a 6D subspace described by 
{(r,r',f 3 (r,r'), f w (r, r')) : {r, r'} € P}, which cor- 
responds to the physical state in which the position of 
two electrons fixes the positions of all the others. How- 
ever, we have to keep in mind that we can affect M only 
by varying the one-body potential v(r). Physically, this 
means that t>(r) must compensate the repulsive forces of 
the other N — 1 electrons acting on the electron in r. If 
the position of the other N — 1 electrons only depends on 
r [as in Eq. (18)], then it may be possible to find such a 
v(r). If, instead, we require the degeneracy of M to be 
higher, finding such a v(r) seems a daunting task. 

We shall go through these physical considerations more 
in detail in Sec. IV, where we analyze the solution (18) 
and write down an explicit equation for v(r). Before do- 
ing so, we first show that we can consistently construct 
an antisymmetric wavefunction of the form (13) such 
that |"0(ri, r 2 )| 2 is a distribution that is zero every- 
where except on the minimizing set M. 



III. ANTISYMMETRY OF THE 
WAVEFUNCTION 

We want to construct an antisymmetric wavefunction 
of the form (13) such that its spatial part |V'( r ii ■ ■■,?n)\ 2 
is a distribution that is zero everywhere except in the 
set M where E pot has its global minimum. Given a 
configuration (ri,...,rjv) of M, we first notice that all 
its electronic positions must be different from each 
other: without the kinetic energy operator T, two or 
more electrons on-top of each other would make E pot 
infinite because of the singular nature of Coulomb repul- 
sion at short distances. Moreover, for any configuration 
€ M, also all its TV! permutations are £ M, since E pot 
is invariant under permutation of its variables rj. Then, 
given the value of the wavefunction for one configuration 
(ri, r^v) <5 M, with spins (cti, itjv), its value for the 
permutations tt of this configuration must be simply de- 
fined as 

^(^(l), I"7r(7V))x( fJ 7r(l), ■■■,0'ir(JV)) = 

(-iri>(r 1 ,...,r N )x(a 1 ,...,a N ) (19) 

where ( — l) 71 " denotes the sign of the permutation n. Be- 
cause the Vi in a given configuration are all different from 
each other, and ip is zero everywhere except on M, the 
N electrons are always localized in different regions of 
space that have exactly zero overlap. This means that 
when we consider the square of the spatial wavefunction 
\tp(ri, ...,r7v)| 2 all the cross terms are exactly zero, and 
thus l^l 2 is completely independent of the choice of the 
spin eigenstate x- 

Similarly, for repulsive bosons (with any spin), ^ can 
be symmetrized. Both the symmetric and the antisym- 
metric choice give the same spatial \ip(ri, rjv)| 2 , and 
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thus the same expectation for V cc : when Coulomb re- 
pulsion between the particles becomes dominant with re- 
spect to the kinetic energy, particle-particle overlap is 
suppressed, so that the particles no longer know wether 
they are fermions or bosons. 

The same construction (19) applies to the SCE state in 
which M has the form (18), since all the considerations 
made so far are obviously still valid for any configuration 
in M characterized by a given reP. 



IV. CO-MOTION FUNCTIONS 

We have seen that the potential v(r) of Eq. (9) for the 
a — > oo limit of DFT must be such that the minimum 
of the 3-ZV-dimensional function E pot of Eq. (14) is de- 
generate over the 3-dimensional subspace M of Eq. (18). 
We shall now analyze how to construct such a»(r), and 
thus which properties the functions f n (r) must satisfy, 
and how they can be determined from the density p(jc). 
Since the solution (18) corresponds to a state in which 
the position of one electron dictates the postions of all the 
others, we call the functions f n (r) co-motion functions. 



A. Properties of the co-motion functions 

If we assume that the v(r) we are looking for is a 
smooth potential with a continuous gradient Vti(r), then 
Epot{[v]\ r i, ■■•, r A r ) is a continuous and smooth function 
of its variables r„, except for configurations with two or 
more electrons on top of each other, which, as explained 
in Sec. Ill, cannot belong to the minimizing set M. In 
this case, a minimizing configuration (ri,...,rjv) G M 
must satisfy the stationarity conditions for the function 
£ P ot(M;ri, ...,r N ), i.e., 

' V«(r 1 ) = E£n^r 
< Vt,(r 2 )=££ 2l £=ftr (20) 



If we insert the solution (18) into Eqs. (20) we obtain 

Vu(rj - | r _ fi ( r )|S 

< Vt,(x)i x=f2(r) = E$ 2 |£ff-twi» (21) 



where we have defined fi(r) = r. Now, suppose that we 
have found N— 1 functions fj(r) and a potential v(r) such 
that the first of the Eqs. (21) is satisfied for all r <G R 3 . If 
we now evaluate this first equation for r = £ n (s) and then 
put s = r, we see that its left-hand-side coincides with 
the left-hand-side of the n th equation. It is then easy to 
verify that, if we want also the right-hand-sides to be the 
same, the functions fj(r) must satisfy the transformation 



properties 

{f 2 (f„ (r )),..., f N (f„ (r ) ) } - {fx (r ),..., f N (r) }\ {f„ (r ) } , 

(22) 

i.e., the set of N — 1 functions f»(x) with i = 2, ...,N, 
when evaluated in x = f ra (r) must yield any permuta- 
tion of the set of N — 1 fj(r) with i = 1, N and i ^ n 
(thus including fi (r) = r) . This means that applying one 
co-motion function ffc to one position r„ of a given SCE 
configuration CeM must always yield a position of the 
same C again. If the functions fi satisfy this property, 
then the fulfillment of the first of the Eqs. (21) automat- 
ically implies the fufillmcnt of the other N — 1 equations, 
and thus the stationarity of the solution (18) for any r. 
One has then to verify that such a stationary solution 
is the global minimum of the 3./V-dimensional function 
E pot ([v],ri, ...,r N ). 

B. The SCE external potential 

Once we have found some functions fj(r) that satisfy 
the properties (22), Eqs. (21) provide N equivalent equa- 
tions for the potential v(r), 

™-±^- 

Equation (23) has the clear physical meaning anticipated 
in Sec. II: the potential v(r) must compensate the net 
force acting on the electron in r, resulting from the re- 
pulsion of the other N — 1 electrons at positions fj(r). 

We start now to see that the a — > oo limit of DFT is 
entirely characterized by the co-motion functions fj(r), 
which also give, as a byproduct via Eq. (23), the external 
potential v(r). In the next Subsec. IV C we shall see how 
to determine these functions from the density p(r). 

C. Co- motion functions for a given density 

The SCE problem for a given density p(r) reduces then 
to the construction of the appropriate co-motion func- 
tions f n (r). To do this, we simply use the quantum me- 
chanical meaning of the electronic density. Since in the 
SCE state the position of the first electron determines 
the positions of all the others, the probability of finding 
the first electron in the volume element dr around the 
position r must be the same of finding the n th electron 
in the volume element df n (r) around the position f„(r). 
This means that all the co-motion functions f„(r) must 
satisfy the differential equation 

p(f„(r))df„(r) =p(v)dv, n = 2,...,N. (24) 

In order to construct the co-motion functions we thus 
have to find the initial conditions for the integration 
of (24) that (i) satisfy the properties (22), (ii) yield a 
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smooth potential v(r) via Eq. (23), and (iii) give the min- 
imum expectation of V ee . In Sec. V we solve this problem 
explicitly for spherically symmetric ./V-electron densities. 

Equations (24) can also be proven by explicitly 
constructing the square of the spatial wavefunction 
\tp(ri, rjv)| 2 that is zero everywhere except on M of 
Eq. (18), and by imposing that the expectation value of 

the operator p(r) = J2iLi <>( r ~ r i) 011 l^( r i> — > r iv)| 2 
is p(r). In this way, we also see that the value of 
\ip(r, f 2 (r)..., fjv(r))| 2 in a given SCE configuration e M 
must be set equal to jjp(r). 

D. The value of W x [p] 

After all, maximising with respect to the potential v 
in Eq. (12) seems to be equivalent to constructing the 
co- motion functions f„(r). In terms of these functions, 
the expectation of the operator V ee in the SCE state [see 
also Eqs. (2) and (5)] reads 

N-l n -L M 

where, again, we have used the convention fi(r) = r. 
Equation (25) comes from the fact that, as said in the 
previous Subsec. IV C, for a given r, the configuration 
(r, f 2 (r), f^v(r)) in M has a weight equal to jyp(r), 
and that the only possible electron-electron distances are 

l f i( r ) - fj( r )l> with i ¥"3, and i,j = 1, --^N- 



V. THE SCE SOLUTION FOR N ELECTRONS 
IN A SPHERICAL DENSITY 

As an example, we consider here N electrons in a 
spherical density, and we explictly calculate the functions 
f„(r) and the value Woo[p] for a few atoms. 

Given the symmetry of the problem, we choose a 
spherical-coordinate reference system = (ri,9i,(f>i) in 
which the position of the first electron defines the z- 
axis (#i = 0, <p\ = 0) and the second electron is on 
the xz plane (0 2 = 0). The 2N — 3 relative angles 
{#2, 83, <f>3, On, <Ajv} are globally denoted by ft. The 
total classical energy we want to minimize then reads 

JV 

E P ot([v},n, ...,r N ,Cl) = ^2v(ri) + V ee (n, ...,r N ,fl). 
»=i 

(26) 

Since the external potential does not depend on the rela- 
tive angles fi, we can decouple the stationarity equations 
(20) into an angular part, which only involves V ee , and 
a radial part that also involves v(r). By solving the an- 
gular part for any given ri, ...,rjy we can define a func- 
tion Q(ri, rjv) for the minimizing angles, valid for any 
spherically-symmetric v. These angles are the solution of 
the electrostatic equilibrium problem for N neutral sticks 



of lenghts r\, rjv having the same point-charge q glued 
at one end, and the other end fixed in the origin, in such 
a way that they are free to rotate in the 3D space. 

We now insert the function 0(ri, r^) in the station- 
arity equations for the radial variables ri, 

' v '( r ^ = -^r^ee(ri ) r 2) ...,rj V ,n)| n(rii ___ iT . j(r) 
< v'(r 2 ) = -^V ee (r 1 ,r 2 ,...,r N ,n)\ n{ru rN) (27) 



As explained in the previous Sec. IV, in order to have 
a smooth (e.g., atomic) radial density p(r) we need to 
choose a v(r) that makes these equations fulfilled by a set 
of radial distances (r, / 2 (r), /jv(^)) m which r can take 
any value in the domain in which p(r) =/= 0. These radial 
co-motion functions f n (r) must have the same transfor- 
mation properties of Eq. (22), as it can be verified by 
inspection of Eqs. (27). 

To construct the radial co- motion functions f n {f) for 
a given density p(r) we proceed as follows. Our starting 
point is Eq. (24) that becomes, in the case of spherical 
symmetry, 

4^ f n {r? p(fn(r)) \f' n (r)| dr - 4irr 2 p(r) dr. (28) 

We then define the function A^ e (r) that gives the expected 
number of electrons between and r, 

N e (r) = [ 4TTx 2 p(x)dx. (29) 
Jo 

N e (r) is a monotonously increasing function of r, and its 
inverse A r ~ 1 (y) is well defined for y e [0,7V). By inte- 
grating both sides of Eq. (28) we see that each function 
/„ (r) is of the form 

f N-^N^UiO)) ± N e (r)) < r < ai 

f n {r)={ N- 1 (N e {f n {a 1 ))TN e {a 1 )±N e (rj) 01 < r < a 2 

(30) 

where the upper (lower) sign corresponds to the case 
fn > {f' n < 0), and [0, 01], [ai, a 2 ], [a k -i, a k ] are 
the intervals in which (i) f' n {r) does not change sign and 
(ii) N e (f n ( ai )) T N e ( ai ) ± N e (r) £ [0, N], so that iV^ 1 is 
well defined. 

The SCE problem for spherical densities thus consists 
in finding the signs, the intervals [<ij_i,aj], and the ini- 
tial conditions /„(0), f n (dk) in Eqs. (30) that (i) make 
the f n (r) satisfy the transformation properties (22), (ii) 
yield, via Eq. (23), a v(r) with a continuous first deriva- 
tive, and (iii) give the minimum expectation value of V e e 
when inserted in Eq. (25). 

We shall now see how to proceed in the case of atomic 
densities. We first treat the case of the He atom, and we 
then generalize the solution to TV > 2. 
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A. N = 2 (He atom) 

Two electrons in a smooth spherical density in the 
a — > oo limit have been already considered in Rcf. 7, 
where the SCE state was proposed on the basis of physi- 
cal considerations. Here, we restate the problem in terms 
of the formalism just derived, and we show that the so- 
lution proposed in Ref. 7 was indeed the correct one. 

When N = 2 we have only one relative angle, 2 . By 
solving the stationarity equations for the angular part we 
immediatly obtain 2 = or it. It is easy to verify that 
2 = is a maximum and 2 — it is a minimum for any 
value n and r 2 . By inserting the minimizing angle 2 = ft 
in the stationarity equations for the radial variables, we 
obtain 



v'(n) 

v'{r 2 ) 



(ri+r 2 )- 2 

(n+r 2 )- 2 



(31) 



These equations admit a solution of the kind (r,f(r)) if 
and only if f(r) satisfies the property /(/(r)) = r, in 
agreement with Eq. (22). This means that, among all 
the possible f(r) of the form (30), we can only choose 
the ones for which / _1 = /. One can verify that the 
only possible choices are then 



fir) 
fir) 

fir) 



N~ 1 (N e (r))=r 

N-HN e (f(0))-N e (r)) 
N- 1 (2 + N e (f(0)) 

N-^l + N^r)) 
N-\N e (r)-l) 



r < /(0) 
N e (rj) r>/(0) 

r < JV-^l) 
r>N~\l) 



(32) 
(33) 



(34) 



We call the choice (32) the "breathing" solution: the two 
electrons are always at the same distance from the center, 
opposite to each other. It is a stationary solution valid 
for any density, since f(r) = r satisfies the differential 
equation (28) independently of the choice of p(r). The 
corresponding external potential is v(r) = — The Hes- 
sian matrix shows that this stationary solution must be 
ruled out, since it is a maximum for _E pot ([i;], n, r 2 , 2 )- 
The choice (33) corresponds to a discontinuous v'(r), 
except in the case 7V e (/(0)) = 2, i.e., /(0) = 00. We rule 
out a discontinuous v'(r) because otherwise our whole 
construction is inconsitent: if v'(r) is not continuous it is 
not necesessary for a minimizing configuration of E pot to 
satisfy the stationarity equations (20). As a double check, 
we also verified, in the case of the He atom, that the 
choice N e (f(0)) = 2 in (33) is indeed the one that yields 
the lowest expectation of V ee , also when compared to 
the form (34) that, again, corresponds to a discontinous 
v'(r). The SCE solution for N = 2 electrons in a smooth 
spherical density is thus the one proposed in Ref. 7, i.e., 



f(r)=N-\2-N e (r)). 



(35) 



It is easy to verify that this stationary solution, with its 
corresponding potential, 



v(r) 



dx 



is the true minimum of E pot ( [v] , r\ , r 2 , #2) • 

Besides the mathematical arguments, the choice (35) 
is the most "physical" one: it makes the two electrons 
always be in two different spherical shells, each one con- 
taining, on average in the quantum mechanical problem, 
one electron [7]. This is exactly what we expect for two 
electrons that repel each other infinitely strong, but that 
have to fulfill the constraint of yielding a given smooth, 
atomic- like density p(r). 



B. N = 3 (Li atom) 

When N = 3 we have 3 angles, 2 ,0s,<j>3, but we im- 
mediately obtain 03 = 0: the three electrons must be on 
the same plane, containing the nucleus, to achieve com- 
pensation of the forces (see Sec. IV B). We then find nu- 
merically, for any given ri,r 2 ,r 3 , the minimizing angles 
2 and #3. 

To construct the radial co-motion functions fair) and 
/3(f) we should, in principle, try out all the possible f n (r) 
of the form (30) that satisfy the properties (22), and se- 
lect the ones that yield a continous v'(r) and the lowest 
expectation for V cc . Instead, we use the physical idea that 
was behind the solution for two electrons of Eq. (35): we 
expect that the correct f 2 and f 3 are the ones that make 
the three electrons always occupy three different spher- 
ical shells, each containing, on average in the quantum 
mechanical problem, one electron. Such radial co-motion 
functions read 



Mr) 
hir) 



N- 1 (2-N e (r)) r<a 2 

N- 1 (N e {r)-2) r>a 2 

N- 1 {2 + N e {r)) r< ai 

N-\A- N e (r)) r> ai 



(37) 
(38) 



(1 + ./W 



(36) 



where ai = A^l) and a 2 = N t T 1 (2). They are dis- 
played, for the Li atom using an accurate fully correlated 
density [13], in Fig. 1. 

To better grasp the physics behind such solution, we 
show, in Fig. 2, the positions of the three electrons in the 
SCE state for the same accurate Li atom density [13]. 
The two big circles are the radii of the spheres contain- 
ing, on average in the quantum mechanical problem, one 
electron (ai) and two electrons (a 2 ). The position r of 
the first electron is varied along the vertical z axis a) 
from to ai, b) from a\ to a 2 , and c) from a 2 to 00. 
The other two electrons have distances from the center 
given, respectively, by fair) and fzir), and angular po- 
sitions given by the minimizing angles 2 (r, f 2 ir), fs{r)) 
and 6*3(r, / 2 (r), / 3 (r)). In each panel, the starting posi- 
tion of the three electrons is represented by a full circle 
(•) and their final position by an empty circle (o). The 
dashed curves represent the "trajectories" of the three 
electrons, and the arrows their direction. Thus, for each 
position r on the z axis of the first electron, the po- 
sitions of the other two electrons are completely fixed, 
apart from the permutation symmetry of electron 2 with 
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Li atom 



-fafr) 




FIG. 1: The radial co-motion functions fzir) and /s(r) for 
the Li atom density. Hartree atomic units are used. 



electron 3. We clearly see from this figure, that in any 
SCE configuration the space is divided in three spheri- 
cal shells, which never contain more than one electron. 
These three spherical shells are the same that, in the 
quantum mechanical problem, contain on average one 
electron: the SCE state makes them contain always ex- 
actly one electron, thus suppressing any accidental clus- 
tering ("correlation suppresses fluctuations" [14, 15]). 

Because of the properties of the co-motion functions 
(which are a consequence of the fact that the electrons are 
indistinguishable), the three panels are equivalent. Our 
physical problem is in fact invariant for a rigid rotation 
of any of the configurations corresponding to a given r. 
This means that in panel b) we could, for each value 
of r, rigidly rotate the position of the three electrons in 
such a way that the electron in the inner shell moves on 
a straight line. In this way, we would reobtain a rigid 
rotation of panel a). The same can be done with panel 
c) . Thus, to calculate the expectation of V ee we only need 
to consider the case < r < a\. 

The solution of Eqs. (37)- (38) is thus very reasonable 
from the physical point of view. Moreover, it can be 
easily verified that it satisfies the properties required by 
stationarity of Eq. (22), and that it corresponds to a 
smooth external potential v(r), which can be computed 
via 

= -^ y -( ri ' r2 ' r3 '^'^)lr,/ 2 (r),/3 W , 02W , e 3(r)' 

(39) 

where 0i(r) is a shortened notation for $i(r, /2(f), /3 (?")). 
The potential obtained by integrating Eq. (39) (again 
using the accurate density for the Li atom of Ref. [13]) 
is reported in Fig. 3. Notice that v(r) is finite, with a 
zero first derivative, at r = 0: the electron-nucleus cusp 
in the density p(r) is compensated by the kinetic energy 
for any value of a (including the very limit a — > oo), as 
discussed in Ref. 16. For large r, we have v(r — > oo) = 
— E^l^ a necessary condition to include as a stationary 
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FIG. 2: The positions of the 3 electrons in the SCE state for 
the Li atom density. The two big circles are the projections of 
the spheres containing, on average in the quantum mechanical 
problem, 1 electron (radius ai) and 2 electrons (radius 02). 
The position r of the first electron is varied along the vertical 
z axis a) from to 01, b) from ai to 02, and c) from 02 to 00. 
The other two electrons have distances from the center given, 
respectively, by /2(r) and fs(r) of Eqs. (37)- (38), and angu- 
lar positions given by the minimizing angles 02{r, f2(r), /3(f)) 
and #3(7-, f2(r), fz (r)). In each panel, the starting position of 
the 3 electrons is represented by a full circle (•) and their final 
position by an empty circle (o). The dashed curves represent 
the "trajectories" of the 3 electrons, and the arrows their di- 
rection. The three panels are actually equivalent (see text). 
Distances are in atomic units. 




8 








-2.5 1 1 1 1 1 1 1 

2 4 6 8 10 

r 



FIG. 3: The external potential v(r) of Eq. (9) for the SCE 
state in the case of the Li atom. The asymptotic expansion for 
large r, v(r — ► oo) = — ^F^, is also reported. All quantities 
are in Hartree atomic units. 



state in a smooth density have one configuration in which 
they are both at the same distance from the nucleus, but 
for three electrons this does not happen. One can also 
easiliy check that, again, the "breathing" solution, with 
its potential v(r) = 1S n °t a minimum for the cor- 

responding £; pot ( [v] , n , r 2 , r 3 , 2 , #3 , fa ) • 

We can at this point clearly see the difference between 
the SCE state for a smooth atomic density (Fig. 2) and 
the more familiar Wigner-crystal-like state (for a case 
with three electrons see, e.g., Ref. 17), in which the three 
electrons are localized at the vertices of a triangle at a 
certain distance from the nucleus, say ro. In this latter 
case, the density becomes very peaked around ro, loosing 
any resemblence with an atomic density, and becoming 
more and more similar to Eq. (17). 

C. N = 4 (Be atom) 



configuration the one in which one of the N electrons is 
at infinity The large-r asympotic expansion — 2/r is also 
shown in the same Fig. 3. 

With this potential v(r) we have also verified that the 
stationary solution of Eqs. (37)-(38) is the absolute min- 
imum of E pot ([v], ri, r 2 , r 3 , 9 2 , #3, 4> 3 ). All these checks 
make us believe that Eqs. (37)- (38) are the correct SCE 
solution for the Li atom: it is difficult to imagine a differ- 
ent solution that yields a lower value for the expectation 
of V cc , since our solution is the one that makes the 3 
electrons be always as far as possible from each other, 
without violating the constraint of yielding the smooth 
density p(r). 

We can thus finally obtain the value VFxJp] = 
—2.6030 Hartree for the Li atom, via 

W oc [ P }+U[p} = 

/ 1 dr4irr 2 p(r)V ee (r,f 2 (r),fs(r),6 2 (r),6 3 (r)), (40) 
Jo 

where we have used the fact that, because the three elec- 
trons are indistinguishable, integrating from to 00 is 
the same as integrating three times from to 01. 

As a final remark, we comment briefly on the fact 
that one never obtains the configuration in which the 
three electrons are at the same distance from the nu- 
cleus, at the vertices of an equilateral triangle. The rea- 
son is that the only f 2 (r) and f 3 (r) compatible with 
this configuration (even for only one single value of r, 
say ro) that also satisfy the required properties (22) 
are the ones corresponding to the "breathing" solution, 
f 2 ( r ) = /3M = r , as it can be easily verified by inte- 
grating both sides of Eq. (28) with the initial condition 
r = tq, f n (ro) = fQ. This gives the two possibilites 
f A (r) = r or f B (r) = N- 1 (2N e (r ) - N e {r)): any choice 
for f 2 (r) and fz(r) which contains the function /s(r) 
does not satisfy the properties (22). This is a difference 
with the N = 2 case in which, instead, the solution is 
f(r) = /s(r) with N e (r ) = 1: two electrons in the SCE 



We have now 5 relative angles to consider. As in the 
N = 3 case, we compute the minimizing angular function 
Sl(n , r 2 , r 3 , r±) numerically. 

We then construct the radial co-motion functions, 
hif), H{t), f&ir), following the same ideas used for the 
case N ~ 2 and 3: we divide the space in four spher- 
ical shells, each containing, on average in the quantum 
mechanical problem, one electron. The radial co-motion 
functions that make the four electrons always be in four 
distinct shells arc then 



h{r) 

Mr) 
U(r) 



N- 1 (2-N e (r)) r<a 2 

N- l {N e {r)-2) r>a 2 

/ N- 1 (2 + N e (r)) r<a 2 

\ N-\6-N e (r)) r>a 2 



--N- 1 

1 " e 



(41) 

(42) 
(43) 



where, again, (Zj = A r ( T 1 (z). These functions satisfy the 
transformation properties (22), and are reported in Fig. 4 
for the case of an accurate correlated density [18] of the 
Be atom. 

We can then obtain, using the analogue of Eq. (40) for 
N = 4, the value W^p] = -4.0212 Hartree for the Be 
atom. 

To check our result, we also computed the correspond- 
ing external potential v (r), which, as shown in Fig. 5, is 
continuous and smooth. We also verified that our sta- 
tionary solution is the absolute minimum for the corre- 
sponding Spot ([v] , ri , r 2 , r 3 , r 4 , SI). 

Interestingly, we found that, for the case of the Be 
atomic density, the four electrons always lie on the same 
plane, containing the nucleus. This effect seems to be due 
to the shell structure of the atomic density: by repeating 
the calculation using a simple exponential density, p(r) = 
/3 3 e~P r /2n, we found that the four electrons are often 
non-complanar. Instead, by using other test densities 
with a shell structure, we found again the four electrons 
always on the same plane. The shell structure seems to 
force the four electrons to have distances from the nucleus 
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FIG. 4: The radial co-motion functions fair), fe{r) and /4(f) 
for the Be atom density. Hartree atomic units are used. 



FIG. 6: The radial co-motion functions for the sphericalized 
B atom density. Hartree atomic units are used. 
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FIG. 5: The external potential v(r) of Eq. (9) for the SCE 
state in the case of the Be atom. The asymptotic expansion 
for large r, v(r — ► 00) = — ^jr-, is also reported. All quanti- 
ties are in Hartree atomic units. 



that are distributed in a highly nonuniform way, so that 
forces compensation (see Sec. IV B) can be achieved only 
if they are on the same plane. 



D. The general iV-electron solution 

This way of constructing the solution can be general- 
ized to any number of electrons N. The radial co-motion 
functions that make the N electrons always occupy N 
different shells containing, on average in the quantum 
mechanical problem, one electron, are as follows. Define 
an integer index k running for odd N from 1 to (N — 1)/2, 
and for even N from 1 to (N - 2)/2. Then 



hk+\{r) 



f N-\2k-N e {r)) r<a 2k 

J 2 ^ r >-\N- 1 (N e (r)-2k) r>a 2k 

N- 1 (N e {r) + 2k) r<a N _ 2k 
N~\2N -2k-N e (r)) r > a w _ 2fe , 



(44) 
(45) 
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1 - x_ f 4 (r),X/ f6(r) / > 
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FIG. 7: The radial co-motion functions for the sphericalized 
C atom density. Hartree atomic units are used. 



where aj = N~ (i). For odd N, these equations give all 
the needed N — 1 radial co- motion functions, while for 
even N we have to add the last function, 



f N (r)=N- 1 (N-N e (r)). 



(46) 



Using these radial co-motion functions, which satisfy 
the properties (22) for any N, we calculated the SCE 
value Woo [p] for accurate sphericalized densities of the B 
and C atoms [19], and for the Ne atom density [18, 20]. 
The resulting radial co-motion functions are displayed in 
Figs. 6-8. The results for Woo[p] are reported in Table I, 
and the corresponding potentials v(r) are shown in Fig. 9. 
In all cases we find the TV electrons non-complanar. 



VI. COMPARISON WITH THE PC MODEL 

In the absence of spherical symmetry there seems to be 
no general strategy for finding the initial conditions to in- 
tegrate Eqs. (24), although we are still investigating this 
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FIG. 8: The radial co-motion functions for the Ne atom den- 
sity. Hartree atomic units are used. 
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FIG. 9: The external potential v(r) of Eq. (9) for the SCE 
state in the case of the Ne atom and of the sphericalized B 
and C atoms. All quantities are in Hartree atomic units. 



problem. So far, such cases can only be treated approx- 
imately, using instead of the exact expression (25) the 
PC model (6) for the functional Wx,[/>]. The PC model 
has been tested successfully [8] against the meta-GGA 
functional of Ref. 21 that is expected to be accurate in 
the limit a — > oo. The present exact (SCE) results, in 
contrast, provide a rigorous test for the PC model. 

Table I documents reasonable values from the PC 
model for all the atoms considered here. A severe test 
is the trivial case of the hydrogen atom H, with the one- 
electron density p(r) = ^e~ 2r , where (V^ ) = for any 
value of a and thus Woo [p] = — yg Hartree. For the other 
atoms in Table I, the PC error is less than 63 mH w 39 
kcal/mole. 

The quantity Woo [p] refers to the limit of infinitely 
strong repulsion between electrons. Being the sub- 
ject of the present paper, this extreme limit is not di- 
rectly relevant for real electron systems such as atoms 
or molecules. Rather, it is studied here in order to be 
combined with additional information from the opposite 





W^[p] (H) W™\p] (H) 


error (mH) 


H 


-0.3125 


-0.3128 


0.3 


He 


-1.500 


-1.463 


37 


Li 


-2.6030 


-2.5559 


47 


Be 


-4.0212 


-3.9608 


60 


B 


-5.7063 


-5.6502 


56 


C 


-7.7817 


-7.7192 


63 


Ne 


-20.035 


-20.000 


35 



TABLE I: Comparison of the values Wxjp] in Hartree atomic 
units obtained with the SCE construction, and with the PC 
model [8]. The absolute errors of the PC model are also re- 
ported. 



weak- interaction limit a — > 0. In this way, properties 
of real atoms and molecules, including those with weak 
correlations, can be predicted accurately [6] without deal- 
ing explicitly with the complicated wave function at the 
realistic interaction strength a = 1. Both the extreme 
limits a — > and a — > oo are mathematically much sim- 
pler than the realistic situation at a = 1. Consequently, 
more relevant than the error of the quantity Wx>[p] it- 
self is the error it causes in the ISI correlation energy [6] 
El sl [p] = Elf[p] - E x [p}. Since the ISI model of Ref. 6 
uses, along with W>o[p], also the coefficient W^Jp] which 
is not known exactly, we consider here the earlier ISI 
version SPL of Ref. 5, 



with 



Q[p] 



y/l + 2Q[p] a 
2\E?™[p}\ 



Analytical integration of f Q daW^ PL [p] — E x [p] yields 



(48) 



Ef L [p) = (E x [p}-W 00 [p}) 



y/l + 2Q[ P ] - 1 

Q[p] 



- 1 



(49) 

The correlation energies E^ PL [p] , evaluated with the val- 
ues Wqo [p] of Table I, are reported in the columns "(PC)" 
and "(SCE)" of Table II. While the SPL prediction 
comes quite close to the exact correlation energy of the 
He atom, the case of the Be and Ne atoms is less satisfy- 
ing. Nevertheless, the improvement beyond the second- 
order correlation energy is remarkable, since the only ad- 
ditional information used here is the coefficient Wxjp] 
(but not W^[p}). 

For the He, Be, and Ne atoms, the error introduced 
by the PC model in E^ ph [p] is always less than 1 mH, 
showing that in the case of neutral atoms the SPL-ISI 
correlation functional is not too sensitive to the exact 
value of Woo [/:>]• However, we have to keep in mind that 
neutral atoms are systems that resemble much more to 
the noninteracting KS system (a = 0) than to the SCE 
state (a = oo). We expect that for more correlated sys- 
tems (e.g., streched bonds) the ISI correlation energy is 
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EM 




(PC) 


(SCE) 


E« xact [ p ] 


He 


-1.0246 


-0.0503 


-0.0413 


-0.0418 


-0.042 


Be 


-2.674 


-0.125 


-0.1054 


-0.1061 


-0.096 


Ne 


-12.084 


-0.469 


-0.4205 


-0.4207 


-0.394 



TABLE II: The exchange energy -E x [p], the second order 
Gorling-Levy correlation energy £ c GL2 [p] [22], and the esti- 
mate of the correlation energy E c [p] from the SPL ISI model 
[5] using the values W^ E [p] (SCE) and W* c [p] (PC), com- 
pared with "exact" values. All energies are in Hartree atomic 
units. 



much more sensitive to the exact value of [p] . The 
investigation of such cases will be the object of future 
work. 



VII. CONCLUSIONS AND PERSPECTIVES 

We have reformulated the strong-interaction limit of 
density functional theory in terms of a classical problem 
with a degenerate minimum, obtaining a consistent solu- 
tion for this limit. Even if an antisymmetric wavcfunc- 
tion can be explicitly constructed, the strong-interaction 
limit of DFT is entirely characterized by much simpler 
quantities, the so called "co-motion" functions, which are 
related to the electronic density p(r) via the differential 
equation (24). The results of this work can be useful 
to construct interpolations between the weak- and the 
strong-interaction limits of DFT, and to test other ap- 
proximate functionals in the a — > oo limit [23]. 

Future work will be mainly devoted to the calculation 
of this limit for non-spherical densities, and to study the 
generalization to the next leading term, W'^p]. The cal- 
culation of the intracule density (the probability distribu- 
tion for the electron-electron distance) will be also carried 
out in connection with the correlation energy functional 
constructed in Ref. 24. 
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APPENDIX A: HARMONIC INTERACTIONS 

It is instructive to see how the SCE solution changes 
if we replace the Coulomb repulsion with an harmonic 
interaction. Since in this case we have the exact solution, 
we can also clarify how the SCE construction presented 



in this paper gives the a — > oo limit of DFT. We thus 
consider the exactly solvable 2-electron hamiltonian [25] 



4(V? + V^) + ^H+r 2 2 )-||r 1 -r 2 | 2 , (Al) 



and we analyze the case in which a — > oo and the density 
is kept fixed, equal to the one at a = 1, by means of a 
suitable k a . 

The^ ground-state one-electron density of the hamilto- 
nian H a is given by 



2^k a (k a -2a) 
\/k a — 2a + 



(A2) 



so that k a must keep f3 independent of a and thus satisfies 
the equation 



y/k a (k a - 2a) y/hih-2) 



\/k a — 2a 



Vfci - 2 + 



(A3) 



where k\ > 2 in order to have a bound system. We easily 
find 



ff 2 



k a ^ 00 = 2a+ , — + 0{a- 1 l 2 ). 



(A4) 



The square of the spatial wavefunction is equal to 

3/2 

|^a(ri,r 2 )| 2 = 



y/k a (k a - 2a) 



e 2 



(A5) 



Now, we want to compare the exact result with the 
SCE construction. We thus follow the same steps done 
for the case of the He atom in Subsec. V A. Again, we find 
that the minimizing angle is 9 = n and that the radial co- 
motion function must satisfy the property f(f(r)) = r, 
so that the only possible choices are the ones of Eqs. (32)- 
(34). In this case, however, the "breathing" solution of 
Eq. (32), f(r) = r, is the one that yields the minimum 
expectation of U V CC " , and that is the global minimum of 
the corresponding E pot . The very different nature of har- 
monic forces with respect to Coulomb repulsion makes 
the solution completely different also from a qualitative 
point of view. The "breathing" solution, that was a max- 
imum for Coulomb forces, becomes the true minimum in 
the case of harmonic forces. In this latter case, in fact, 
the interaction is not singular at zero electron-electron 
distance so that the two electrons can get on top of each 
other (and they do, at the "nucleus"). This means that 
in this case even in the a — > oo limit the spin state can 
still play a role. The SCE one-body potential is obtained 

by 



v'(r) = r + f(r) = 2r, 



(A6) 



so that v(r) = r 2 agrees with the leading term (0(a)) 
of Eq. (A4) (remember that in the SCE construction we 
directly consider all the quantities divided by a). 
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In order to compare the SCE solution with the ex- 
act wavefunction of Eq. (A5) in the a —* oo limit, we 
have first to keep in mind that Eq. (A5) is the solution 
of the hamiltonian H a , so that it includes also the ki- 
netic energy operator T. Equation (A4) shows that the 
external potential V^t that keeps the density fixed has 
the large-a expansion V£ t = aVoo + V + 0(a -1 / 2 ). The 
SCE construction only gives , but the exact solution of 
Eq. (A5) always includes also the next term, Vo, because 
it is of the same order of T. The comparison with the 
SCE solution can be done in the following way. Replace 
in Eq. (A5) the large-a expansion of k a up to orders a , 
k a — ► 2a + We can rearrange the result, with c = 



When a — > oo, the first gaussian in Eq. (A7) with its 
factor in front tends to the distribution <5(ri + r 2 ). The 
second gaussian gives ip(ri). Thus, as a — > oo we have 



\ip a (r 1 ,r 2 )\' 2 



\/2a + C- i ( y2H+H-^)|r 1+ r 2 | 2 
2lT ) 

3 / 2 _ which is exactly the SCE solution, since in the case of 

e _v ^ (r i +r 2( ) A7harmonic forces f(r) = -r. 



lim |i/> a (ri,r 2 ) 

a.— >oo 



5(n+r 2 ), (A8) 



1 / 2^c(2a + c) 



7T 3 / 2 I y/2a + C - Vc 
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